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Abstract 

Iti this report, we present comparisons of results of direct numerical 
simulations of turbulence with both laboratory data and self -similarity 
theory for the case of the turbulent wakes of towed, axisymmetric bodies. 
In general, the agreement of the simulation results with both the labora- 
tory data and the self-similarity theory is good, although the compari- 
sons are hampered by inadequate procedures for initializing the numerical 
simulations . 
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1. Introduction 

This paper reports the use of direct numerical simulations to study 
the wakes of towed, axisymmetric bodies. Most past numerical work on 
turbulence has employed the Reynolds-averaged equations of motion, which 
are necessarily supplemented by a series of ad hoc assumptions. However, 
direct numerical simulations of turbulence involve the numerical solu- 
tion of the fundamental equations of motion to obtain the detailed 
structure of the turbulent field. What is obtained from the calcula- 
tions, then, is the time history of the very complex turbulent velocity 
field. Tills method is analogous to performing experiments in the 
laboratory, and statistical results can be obtained by temporal, spatial, 
and/or ensemble averaging over the computed flow field. 

Some of the earliest work using direct numerical simulations of 
turbulence was applied to low enough Reynolds numbers (R-^ <50, see, 
e.g., Orszag and Patterson, 1972) so that all the dynamically significant 
scales of motion were adequately resolved by the numerical algorithm. 

For higher Reynolds number flows, however, there is no hope to resolve 
all the relevant scales of motion using computers available either now 
or in the near future. Thus, for high Reynolds number flows, the meteor- 
logists’ approach to numerical weather prediction is followed. In this 
approach, the smaller, subgrid-scale motions, which cannot he resolved 
on the computational mesh, are filtered out of the fundamental equations. 
Only the grid-scale motions are computed directly. (For this reason, 
when direct numerical simulations are applied to higher Reynolds number 
flows, the technique is sometimes referred to as large eddy simulations.) 
This filtering process introduces subgrid-scale effects in the form of 
stresses which are analogous to Reynolds stress terms appearing in the 
Reynolds-averaged equations of motion. 

The subgrid-scale stresses are normally modeled in a manner anal- 
agous to the modeling of the Reynolds stresses. There is reason to 
believe, however, that the subgrid-scale modeling will have a better 
chance of success than conventional turbulence modeling. First, the 
smaller scale turbulent motions behave in a more universal manner than 
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tbe larger scale motions. Hence, it is more likely that a single model 
for the subgrid-scale stresses may be applicable to a broad class of 
flows, Second, the subgrid-scale motions usually contain only a small 
portion of the turbulence energy, so that the modeled terms have less 
influence on the overall dynamics than the Reynolds stresses. Thus, 
this approach may overcome some of the drawbacks to the conventional 
approach to computing turbulent flows using the Reynolds-averaged 
equations. This technique, however, will always possess the disad- 
vantage that a more difficult numerical problem must be treated, which 
is necessarily three-dimensional and time-dependent, and for which 
greater spatial and temporal resolution is necessary. 

As mentioned, this procedure is analogous to performing experiments 
in the laboratory. However, it has the advantages that (i) much more 
statistical information of interest can be obtained (since the entire 
flow field is known at every time step) , (ii) parameters can be varied 
easily, and (iii) experimental conditions are more controllable. The 
technique also offers the advantage of circumventing the closure prob- 
lems that have plagued turbulence theory. Because of these advantages, 
direct numerical simulations of turbulence can serve a variety of pur- 
poses. For example, they can be usefully employed to test theories from 
the more exotic (e.g., the direct interaction approximation) to Reynolds 
stress closure models. They can also be very useful in obtaining a better 
understanding of the underlying physical processes. And these simula- 
tions may be ultimately used as a predictive technique for applied 
problems. However, before this approach can be used with confidence, 
especially when applied to high Reynolds number flows, its validity and 
accuracy under various conditions need to be determined through careful 
comparison of the results of the simulations with experimental data. 

Some work has been carried out with the purpose of establishing the 
validity and accuracy of numerical simulations of turbulence. For 
example, research at Stanford and NASA Ames ■ (Mans our at al. , 1977; 

Fersiger et al., 1977) has led to the successful prediction of the 
large-scale features of homogeneous turbulence decay. There has also 
been some success in the comparison of simulation results with channel 
flow data (Deardorff, 1970; Grot 2 bach and Schumann, 1977). However, in 
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these studies the flow in the region near the wall, where many of the 
physical processes of current interest occur, was modeled. Recent work 
by Moin et al. (1978) has examined removing this latter limitation by 
computing the flow in the entire boundary layer. While significant 
contributions have been made in establishing the validity and accuracy 
of numerical simulations,, it is apparent that much more work needs to be 
done to establish the capabilities of the method. 

The study of free shear flows may offer an excellent vehicle for 
determining the capabilities of direct numerical simulations to treat 
turbulent shear flows. While possessing many of the dynamic processes 
of all turbulent flows, free shear flows are generally easier to treat 
numerically than bounded flows. This is because, first of all, the 
boundary conditions are usually easier to implement numerically, and 
secondly, the resolution requirements are much less stringent. In 
particular, the ;.*.v?erful pseudospectral numerical methods can be easily 
applied to some tree shear flows. 

We have chosen to study the wakes of axisymmetr ic , slender, towed 
bodies. This choice was made for several reasons. First, there is 
extensive data available for this case (e.g., Chevray, 1968; Pao and 
Liu, 1973; Bukreev et al. , 1972). Secondly, as long as (i) the mean 
velocity deficit is small with respect to the body speed, and (ii) the 
lateral scale of nonhomogeneity is small with respect to the axial 
length scale, then it can be shown (see Section 2.3) that this flow can 
be treated as a time-dependent problem (with time interpreted as axial 
variation in the data) which is statistically homogeneous in the axial 
direction. This leads to significant simplifications in the numerical 
procedures, and also gives us the opportunity to study a statistically 
time-dependent flow. The two conditions described above are satisfied 
in the late wake in the experiments. Another advantage of late wakes of 
axisymmetric bodies is that they appear to be approximately self-preserv- 
ing, so that the flow develops in a fairly simple manner. For example. 
Figure 1 is a plot from Pao and Liu (.1973) showing the collapse of 
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several sets of data in the late wake, and also the agreement of the 

J* 

data with the self -preserving decay rate". 

The main objective of the work presented is the determination of 
how well direct numerical simulations can model a simple turbulent shear 
flow. Our method of approach is to initialise our calculations with 
available walce data at a specific location downstream of the body, and 
then compute the subsequent development of the walce. The results of the 
numerical simulations are then compared with data downstream of the 
original plane to determine how well the wake has been simulated. The 
effects of a subgrid-scale model Cor lack of one), of the initialization 
procedure, of the numerical algorithms and other features of the simula- 
tions are examined. 

In the next section we discuss the techniques that we used (e.g., 
numerical methods and initialization procedures) to simulate axisym- 
metric wakes. In Section 3 we present the results of our simulations 
and comparisons with laboratory data. Section 4 contains a discussion 
of alternative initialization procedures. Finally, in Section 5 we 
present our conclusions and recommendations. 


A In this plot, U is the maximum axial mean velocity, 

U o the body speed, x is axial distance measured from the 
body stern, D is maximum body diameter, is the drag 


TJ D 
o 


coefficient, and ■■= ---■ , where V is the kinematic viscosity. 


V 
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2. Method of Approach 

Before examining the results of our simulations, it is important to 
first understand the methods used to perform the simulations. These 
methods involve specifically the numerical techniques, the subgrid-scale 
model, the x -*• t transformation used to relate the simulations to 
laboratory data, and the initialization procedures. These topics will 
he addressed in this section. 


2.1 Numerical Methods 

The numerical methods used to solve the Navier-Stokes equations are 
the same as employed by Orszag and Pao (1974) . The Navier-Stokes equa- 
tions are written in rotational form: 




« - tixW 

<v w 



7f> p + v V (& 



( 1 ) 



where to = Vxu is the vorticity vector, and (3) 



is the dynamic pressure. This form of the equation ensures energy 
conservation in our numerical scheme. 

The initial conditions will, be discussed in detail in Section 2.4. 
The wake is chosen to be situated so that the axial direction, the 
direction of non-zero mean flow, is taken to be the x^ direction (Fig- 
ure 2). As discussed below, the flow will be assumed to be statisti- 
cally homogeneous in this direction. Periodic boundary conditions have 
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been successfully employed in the past along directions of statistical 
homogeneity (see, e.g., Orszag and Patterson, 1972), so we choose 
periodic boundary conditions in this direction. In the radial (x^ and 
x^) direction, there is no mean flow (see Section 2.3) so that inflow- 
outflow conditions are not necessary there. Thus, for convenience, and 
also to employ pseudospectral numerical methods most easily, periodic 
boundary conditions are also used in the lateral direction (x^) and 
periodic no-stress (free-slip) conditions in the vertical (x^) * The 
latter conditions are used because the computer program also was written 
to treat density stratification effects in the vertical direction. 

sV 

Pseudospectral methods are used to solve the governing equations . 

The equations are solved on a 32 x 32 x 33 grid in physical space. The 
pseudo-spectral method involves the evaluation of the derivatives using 
the Fourier expansions, i.e., if 


u(x) = X. u (k) e 

IkK# 


lkx 


then is evaluated as 


&7t 


= y ix u(k) e 

»*ii i * 

IkKsK 


ikx 


(5) 


( 6 ) 


The method is similar to the spectral method except that aliasing effects 
which arise in the computation of the nonlinear terms are not excluded 
(e.g.. Fox and Orszag, 1973). The method appears to be as accurate as 
the spectral method, however, and is roughly a factor of 2 more 
efficient computationally. 

The computations presented herein were performed on the NASA Ames 
CDC 7600. The calculations require about 3 sec/time step for the 32 x 
32 x 33 computational mesh. The numerical algorithm uses fast Fourier 


^Pseudospectral methods are approximately twice as accurate in each 
space dimension as comparable finite different schemes with the same 
resolution, and are competitive in computational efficiency for problems 
with simple geometry. Thus, they offer significant advantages when 
applied to problems of the type described here (Orszag and Israeli, 
1974). 
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transform methods to evaluate transforms and inverse transforms, and 
uses leapfrog time differencing on the nonlinear terms, and Crank- 
Nicolson (implicit) time differencing on the viscous terms. 

The numerical simulation data presented is mostly in the form of 
spatial averages. To perform these averages, we take advantage of the 
axial homogeneity and axisymmetry and divide the flow field into concen- 
tric annuli with axes along the x^ axis (Figure 3) . The gap width of 
each annulus is taken to be the mesh spacing (Ax, ) . The average value 
of a quantity q in the annulus is obtained by summing the values 
of q for each grid point in the annulus, and then dividing by the 
total number of points in the annulus, e.g.. 


Ave.M^£ (y) 



(7) 


where N is the number of grid points in the n — annulus. In Table 1, 
n 1 

we show N for each n , and also - which is roughly proportional 

“ K 

to the statistical error expected in computing an average. Note that 
the error is expected to decrease substantially with increasing distance 
from the axis, until the annuli intersect the boundaries of the 

computational domain. 

2.2 Subgrid-Scale Modeling 

The basic idea of direct numerical simulations of turbulence is to 
numerically compute the detailed evolution of the complex turbulent 
velocity field. Because the simulations are inherently three-dimensional, 
the best possible numerical resolution using present-day computers would 
span about one to two decades in wave number space. However, the laboratory 
data that we choose to model (Pao and Lin, 1973) has significant inter- 
actions occurring over many decades. It is obvious that the simulations 
cannot treat all the details of this flow. 

In order to handle the smaller scale motions which cannot be resolved 
in the computer simulations, subgrid-scale modeling is usually intro- 
duced. This method was first suggested by Smagorinsky et al. (1965) 
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for general circulation models of the atmosphere, while Lilly (1967) 
first studied its use in turbulence simulation. To employ subgrid-scale 
modeling, the Navier Stokes equations, (1) and (2), are first filtered 
to eliminate the smaller scale motions. Generally a spatial average is 
used, of the form 

«(&'*>*) dx' , t8) 

•'R 

where G is a filter function and R encompasses the entire flow 
domain (see, e.g., Leonard, 1974). Here we have used the symbol <u> 
to denote the filtered value of the vector u . The flow variables can 
thus be expressed as a filtered part, or large scale part, <u> , and a 
subgrid-scale component, say u 1 . Thus 


ii tL ( L( ) + tl* . 

When (9) is substituted into equations (1) and (2) , and the filtering 
operation is performed on the equations, there results 

&<«*> + &<<*><“■ i>> = - fi<f> + 

■4-[ <u: i u i> + < <t( ^ a \ > + < u i <u i>> ] 

i.<^> =° 


Note that, in general, «u>> ^ <u> . 

Equations (10) and (11) are the dynamic equations for the resolv- 
able (large eddy) motions. The last term in equation (10) represents 
the effect of the filtered (sub grid-scale) motions on the large scale 
motions. Clearly a closure problem, analogous to that resulting from 
Reynolds averaging of the original equations, has resulted from the 
filtering process. 
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la order to proceed further, the effect of the subgrid-scale motion 
on the resolvable scales needs to be modeled. Smagorinsky et al. (1965), 
realizing the analogy with the Reynolds averaging, modeled the subgrid- 
scale stresses using the eddy viscosity concept. Deardorff (1974a, 

1974b) generalized this by introducing equations for higher order subgrid- 
scale quantities, e.g., for the subgrid-scale kinetic energy and subgrid- 
scale stresses. Leonard (1974) put this work on a firm theoretical 
foundation by carefully defining the averaging procedures and separating 
out the different effects. In more recent years, there has been consider- 
able activity in this problem area, from finding theoretical justification 
for the models (e.g. , Love and Leslie, 1977). to careful comparisons with 
data (e.g., McMillan and Ferziger, 19.78). 

Qrszag (1974) has suggested an alternative to the subgrid-scale 
procedure which might he applicable to free shear flows. It is based on 
the concept of Reynolds number similarity, which is that if the Reynolds 
number of a free shear flow is high enough, then the statistical proper- 
ties Which depend on the large-scale features of that flow (e.g., turbu- 
lent energy, Reynolds stress), will be independent of Reynolds number. 

As the Reynolds number of a free shear flow is increased, the wave 
number range over which viscous dissipation occurs becomes widely sepa- 
rated from the wave numbers containing the bulk of the turbulent energy. 
(Figure 4.) For example, the ratio of the dissipation scale. 



( 12 ) 


(Here u 1 is a characteristic turbulence velocity and & is a charac- 
teristic length scale of the energetic turbulence.) Generally, is 

several orders of magnitude or more so that the energy and dissipation 
scales are widely separate. It is argued that the nonlinear inter- 
actions are somewhat local in wave number space, so that the large 
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eddies are not directly affected by the dissipation scale motions. Thus 
the larger scale motions are essentially inviscid, and transfer their 
energy to smaller scales, to be ultimately dissipated by the much smaller 
scale motion. Thus the large scale behavior of these flows should be 
inviscid, independent of the Reynolds number. Hence, the term Reynolds 
number similarity, 

Reynolds number similarity has been experimentally established for 
the decay of grid turbulence (see, e.g., Townsend, 1976). Also, recent 
experiments on turbulent shear layers has lent support for this principle 
(see Roshko, 1976, for a review of this work). Visualisation studies 
have shown the qualitative similarity of the large-scale structure of 
these flows for a broad range of Reynolds numbers. Quantitative data 
from these experiments also supports the principle. Furthermore, the 
application self-similarity theory to turbulent wakes, which implicitly 
includes the Reynolds-number-independence hypothesis, gives results 
which agree well with laboratory data. This last point will be discussed 
in more detail in section 2.3. 

Past numerical simulations of turbulence which have not used 
subgrid-scale modeling, but which have resolved all the dynamically 
significant scales, have been successful for values of R^ up to about 
50, where R^ is the Reynolds number based upon the Taylor microscale 
\ . (This is still a fairly low-to-moderate value for laboratory studies 
For example, Chevray ' s (1968) Reynolds number was about 700, Pao and 
Lin’s (1973) was about 90, and Champagne, et al. (1979) reported a value 
of 130 in their homogeneous shear flow experiments.) It is possible 
that this Reynolds number is high enough for the flow to be approxi- 
mately Reynolds number independent. We concluded that this was an 
important possibility to address. Thus, in the work presented in this 
report, no subgrid-scale model was used. The R^ was initialised to be 
roughly 50, and we hoped to conclude whether this was large enough for* 
Reynolds number similarity to hold. We plan in the future to recompute 
these cases using a subgrxd-scale model, in order to determine the 
difference in the two approaches. Finally, one might consider our 
approach, i.e., assuming Reynolds number similarity and computing a flow 
with much lower Reynolds number, but with all the scales resolvable, as 
a very simple sub grid-scale model. 
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2 .3 Self-Similarity 

In order to examine axisymmetric wakes further, It is useful to use 
cylindrical polar coordinates- We use the (r,6,x) coordinate system 
shown in Figure 2, where r is the radial coordinate, 0 is the circum- 
ferential coordinate, and x the axial coordinate (measured from the 
body stern), taken to he positive in the direction opposite that of the 
body motion. Let (u , Ug , u ) denote the corresponding velocities. 
Then equations (1) and (2) are 


~r^ u r~^p + a x S< dr - -’p | 

3V ') + - 'p |e u ® * 55?. 




It U S ■» «r|r% + ^ I© ^6 + + U*4. U {> - - 


ftS +• + ^t|b Uk + U Xt&. U fc 

( r Jr 1 *-) + h 5s 1 ■* 


+ l ! ii 


p3x.« > 


1 i frit \ + 4,% + 1 « v 
T r r Q r ax y - 


0 


(13) 


(14) 


(15) 


(16) 


In order to obtain the Reynolds-averaged equation of motion, we 
decompose the velocity vector and pressure into (ensemble) mean and 
fluctuating parts, i.e. , 


(a, ,a s ,0^.) - (yj+us , V+ v ,U+ u.) 


■» P t ft . 


US) 
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Subs tituting (17) and (18) into (13) to (16), ensemble averaging, and 

assuming statistical stationarity and axisymmetry, and no -swirl (V=0) 

* 

gives 

y/ f r W + U *- p irT* -<-vjjr(r s^ x '^) W] (19) 

V irV + V 3x.V [y il* 5? ) - Jjj* V ] CM) 

'7 Tr^ m ) “ ^ 

~2 ~2" 

Equations for the higher order quantities (e.g. , u , uw , w ) can be 

obtained by subtracting equations (19) and (20) from (13) and (15) to 

obtain equations for the fluctuating components, multiplying by w, v, or u 

and ensemble averaging (see, e.g., Tennekes and Lumley, 1972). Statistical 

stationer ity implies that we have fixed our coordinate system to be 

moving with the body, at a speed relative to the fluid which we call 

U . 
o 

We are interested mainly in the far wake of an axisymmetric body 
where the flow is approximately self-similar. In this region the 
following conditions are approximately satisfied (see Tennekes and 
Lumley, 1972, for a more detailed discussion of self-similarity) : 

i. The axial variation of any statistical quantity (q) is much 
smaller than its radial variation, i.e., 

« l-T 


*An overbar ( ) will be used throughout this text to denote an ensemble 
average. 


X35 
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Ay >/' Ap y 

where Z and Z are charact eristic length scales of axial and radial 
ir 

variations in mean quantities; 


ii. The maximum value over which the mean velocity varies, say 
U ■, is much smaller than the body speed, 

TJ £ 

ir> ir 

iii. The Reynolds number of the flow, — “ , is large enough so 

that the viscous diffusion terms in equations (19) and (20) 
can be neglected with respect to turbulence diffusion terms. 

With these conditions, equation (20) reduces to approximately* 

“0. C22) 


Following the self-similarity theory, we assume the axial mean velocity 
and the Reynolds stress to be of the form** 


IftqxV = XSJiA k (n) 


(23) 

(24) 


W ^r 

*Scaling equation. (21) results in — ~ , so that W is small with 

respect to U in the late wake. x 

**ln everything that follows , we take U to be the axial mean velocity 
deficit, i.e. , the difference between the axial mean velocity, measured 
in a coordinate system moving with the body, and the free stream velocity 

U . ,. ■ ■ . . ■ . 

o 
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where f| = is the similarity variable. The velocity scale 

and the length scale £, are to he determined. Mien (23) and (24) are 
used in (22) , and the further condition derived from (22) , 


/ 

V rV Jr = cor si em? 


(25) 


is imposed, there results 


^ Ax 




(26) 


£ M -By 


% 


(27) 


where A and B are constants to be determined empirically. There- 
fore, for the flow in the wake to remain self-similar, the velocity 

-2/3 

scale (e.g., peak velocity) should decay a? x , and the length 

1/3 

scale should increase as x . 

To go any further, a closure assumption is needed. Following 
mixing-lengtli theory we let 


itur =a. - Up ^ *> 

VV- C SUm, 


(28) 

(29) 


and C is a constant to be determined empirically. Then equation (22) 
can be solved to give 


u 


(\) = ax/’ ( - & 

) 


8 

iA n e 


4 AC ' ^ ) 

xt> (- ~£— 

*r\ fAC 




(30) 

(31) 


*Note that this closure assumption does not include dependence of U w 
on viscosity, i.e., Reynolds number similarity has been implicitly assumed. 
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Equations (26) , (27), (30) and. (31) are the main results from the self- 
similarity theory. 

Defining the total crossplane mean kinetic energy E hy 


f m V 

E — J h J ^ Z 


(32) 


o a 

and the total crossplane turbulent energy e by 

rZW & 

e ~ l r< 


r E / 

J& ( dr -|r 


u 2 + v e ■+ 


7 


(33) 


o o' 

then self-similarity theory can also be used to determine the axial 
behavior of E and e . Substituting (26) and (27) into (23), and 
(23) into (32) , and integrating gives 

* . s - ■ : ' "2 A . 

£(x) = a? , 


(34) 


where D is an empirical constant. Following a similar procedure for 
e yields 

(35) 


e(x) - F x 


where F is another empirical constant. 

The approximate validity of the self-similarity theory can be 
established by comparing conclusions from the theory with experimental 
data. Figure 5 is a plot from Pao and Lin (1973) showing the growth of 
the mean velocity wake radius (r ) for several different sets of data 

m ft 

for the wakes of axisymmetric bodies. There is fairly good agreement 
with the theory for each set of data. Figure X is a plot of the decay 
of the maximum mean velocity (U^) for the same data sets (again extracted 
from Pao and Lin) . Again the agreement with the theory is good. Figure 6 
gives the radial behavior of the axial mean velocity in Pao and Lin f s 
data. The solid line is the self-similarity prediction, equation (30) , 
using an eddy viscosity assumption. The agreement of theory and experiment 


*The mean wake radius r is defined here as the distance from the mean 
velocity peak (at r = 0) m to the point at which the . mean velocity attains 
half its peak value. 
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here is also good. Similar results have also been demonstrated (e.g., 

Pao and Lin) for the behavior of the axial component of the turbulence 
~~2 T 

intensity u , the radial component w , and the Reynolds stress 
uw . It is apparent that self-similarity theory is approximately valid 
in the far wakes in these experiments. 

The existence of a self-similar far wake region, and the implied 
conditions necessary, provide the basis for the present calculations. 

In the far wake, assuming conditions (i) and (ii) are satisfied, the 
statistical equations can be obtained from the original equations by 
neglecting all derivatives with respect tc x , except when multiplied 
by IT , Thus the axial mean momentum equation (20) reduces to 


The Reynolds averaged equation for the axial component of the turbulence 
intensity is _ . "TTT 

wi-C?) Kx ) * r s ('' 


— \) 



% X 






mx 

Mb l 


a. 





i. 


With the far-field scaling, this reduces to 


vMi) 


+ 


(, 

±iU(S)) if +/¥?] 

ran TrUi| r a W/ J * 


(37) 


(38) 


Similar equations can be obtained for the other turbulence quantities. 

If is interpreted as time, it is easy to see that (36) and 

o 

(38) are identical to the corresponding equations for a time dependent, 
axisymmetric flow xrLthout mean swirl , which is statistically homogeneous 
in the axial (x) direction. This can be shown to be true for all 
other statistical equations, and is a result of the fact that, in the 
far-field approximation, x derivatives are only retained when multi- 
plied by U q . 
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It is much easier to numerically treat a problem which is statis- 
tically homogeneous in the axial direction and time-dependent, than to 
compute a flow which is statistically steady, but nonhomogeneous in the 
axial direction. In the latter case we have to deal with complex up- 
stream inflow and downstream outflow boundary conditions* Also in the 
far field, the mean radial velocity is negligible with respect to the 
mean axial velocity, so that we need not be concerned with radial inflow 
boundary conditions. Thus, in the present simulations, we choose to 
compute the statistically axially homogeneous, time-dependent problem. 
Time, in our calculations, is interpretable as — in the laboratory 
data. ° 

Note that if the far-field assumptions break down, i.e., if is 

not small with respect to U q , or if the radial scale (J^) is not 
small with respect to the longitudinal scale (Z ) , then it is not 
proper to replace the steady-state wake with a horizontally homogeneous, 
time-dependent wake. This is particularly true in the near field where 

V m ^ cind 

Another way to see the justification for considering the time- 
dependent problem is the following. Fix the coordinate system with 
fluid. Shortly after the body goes by, the mean convection terms (e.g., 




are important. Later on, however, as the mean 
velocity decays, they become negligible. It is then possible to compute 


the time history of a slab of fluid, and, when comparing with data, 

Z 

is small in the late wake, axial 


treat time as 


x 

IT 


Because 


Z 


o s 

derivatives of mean quantities can be neglected. 


2.4 Initialization Procedure 

To initialize the velocity field, we use the data of Pao and Lin 
(1973) at x/D = 36 . This downstream location is in the early stages 
of the approximately self -preserving region. We first separate the 
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velocity field u into a mean component and a turbulence component as 
given by equation (17) , i.e., 

^ Cl7) 
In the self-preserving regime, the only nonzero component or the 
mean velocity is the axial component, U . An excellent approximation 
to the axial mean T U'icity profile in the data is given by the self- 
similar formula (li'.gure 6), 

,2 


where IT is the maximum value of U , and r is defined by 
m m J 


(39) 


We use 


i U *» 

(S0 5 ^e^(-- 
wkate. %( “ Vy, 




% 


Jbr/p - 3&v 
¥d 


(40) 


(41) 


(42) 


(43) 


as initial conditions for our mean velocity field. 

To initialize the turbulent velocity field, we use the method 
suggested by Orszag and Pao (1974) . It is an adaptation to nonhomogeneous 
turbulent flows of the method of Orszag and Patterson (1972), which is 
applicable to homogeneous turbulent flows. We will briefly explain the 
latter firsthand then show how the former method is a modification of 
it. ... . 

In order to ensure incompressibility, Orszag and Patterson write 
the velocity field as the . curl of a stream function vector ip , 
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(Note that this greatly restricts the class of velocity fields considered, 
since in general an arbitrary velocity field can be expressed as the sum 
of an irrotational part (V<jO and a rotational part (Vxip) i see, e.g., 
Batchelor, 1967.) The velocity field is expanded in a discrete Fourier 
series, 

A . . UZ'l 

= L C4S) 

ISKA 

and (44) becomes 

u (k) =. i k x t ' « 6 > 

A 

where ij; is the discrete Fourier transform of $ . The energy spectrum 

A *** ~ 

for ip is determined Uniquely by the energy spectrum for u . The 

^ A -v 

Fourier amplitudes for each k are chosen to be Gaussianly distributed 

•v +i+ 

with zero mean, and the standard deviation chosen to produce the proper 
spectrum for !|; , and hence for u . To ensure homogeneity, the Fourier 
amplitudes for differing wave numbers are chosen to be statistically 
independent. For statistically isotropic flows, the components of a 
particular Fourier mode are taken to be statistically independent but 
identically distributed. Note that, because the initial velocity field 
is comprised of a large sum of statistically independent terms, from the 
Central Limit Theorem (e.g., Feller, 1957) we know that it has a joint 
normal probability distribution. 

The velocity field for a particular realization is set up by first 

A 

numerically selecting each Fourier amplitude ijj from its proper ensemble 

— A 

using a pseudo-random number generator. Then ikx^ is formed to give 

A 'v 

u , and finally the inverse discrete transform is performed to give 
u . 

The modification of Orszag and Pao (1974) also starts by first 

A A 

generating the Fourier amplitude . However, is then transformed 
to physical space to produce t|j , multiplied by a form function, J(i*) 
for the case of an axisymmetric wake, and the product is transformed 
back to Fourier space. The cross product with ik is then taken to 
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ensure incompressibility. The velocity field obtained from this series 
of operations is then 


a 


= v^(T + ) . 


(47) 


We choose the stream function field i|i to be statistically isotropic 
Then its energy spectrum can be expressed as 

3 r /. \ (48) 


£ 


ZK 


E u to) , 


where E^(k) is the energy spectrum for a statistically homogeneous, 
isotropic field obtained by the method of Orszag and Patterson. The 
resulting velocity field defined by equation (47) has the following 
moments; 


ic “ v 


« 2 = v 2 = tt‘ 


7? - u? J z 


s= 

( T z 


+ 



LtV 


= LUST 


^ VtcT - 


Jl 

tut K y 


<3 ■ 

£0 


1 2 C = f £ u (k)c!Ki | c = f 4 (*'i • 


(49) 

(50) 

(51) 

(52) 

(53) 


0" °‘ 

Furthermore, although it is difficult to obtain a general expression for 
the one-dimensional axial energy spectrum, this spectrum at r - 0 
can be shown to be identical to - the one-dimensional spectrum obtained . 
from E^(k) . for the homogeneous, isotropic case. 

In order to initialize the turbulent velocity, we need to select 
the proper energy spectrum. and form function J (r) . We choose the 

energy spectrum E (k) to he given by 


u 


f,.M = £c£A 


a 


k 4 A 4 


(54) 
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Along r = 0 , this produces a one-dimensional axial spectrum for the 
axial velocity u 

a 


Eu (K ) 


l s / ^ 3 Ai A 

W U ° A 




(55) 


and an axial spatial correlation 


aMu(xtt) » 'u d en 


/> (~ s / a ) 


(56) 


which are good approximations to the one-dimensional spectrum and spatial 
correlation in most shear flows. Here A is the longitudinal integral 
scale j i.e.. 




-A = p I afa A$ 3 


(57) 


and u q is defined by equation (53) above. 

From equation (50) and our choice for J (equation (60)), u 


corresponds to the maximum value for u* 1 = ^2 
data this is, approximately. 


V 


In the laboratory 


U 


W 


V, 


a .32., 


(58) 


Yrt 


where u m - max(u') . We choose A to match the corresponding integral 
scale in the data. This is, approximately, 


A 


- v 


(59) 


There is considerable error involved in the estimation of A from the 
laboratory data (as much as 0.3 r ) . Since the proper turbulent length 
scale is of considerable importance in the dynamics of the flow, this 
error could have a significant impact on our results. With the choice 
of u^ and A , our spectrum is completely defined. 
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To choose J'(r) we match our initial radial profile for the axial 
turbulent intensity with the data. A functional form for J which 
gives a fairly close match is 

J~(p) - r /r* ) } (6o) 

. r j 

with — = 2.012. Figure 7 gives the values for the axial turbulent 
m 

intensity u T versus radial distance for two different realizations of 
the ensemble defined above. The solid curve is taken from Pao and Lin 
(1973). Table 1 gives values for a number of relevant statistical 
quantities for one realisation of the initial flow field. 

Although the resulting turbulent velocity field has many of the 
proper statistical characteristics, it is still lacking in many respects. 
First of all, because the Fourier modes were chosen to be statistically 
independent, there is initially not any energy transfer in wave number 
space. This is exemplified by the low values for the velocity derivative 
skewness in Table 1, which are a direct measure of this transfer. It 
has been found, however, that these correlations build up rather rapidly 
(the velocity derivative skewness attains its proper value by roughly 

~~T - 0.5; see e.g., Orszag and Patterson, 1972), and it is not thought 

that this initial lack of correlation is too severe a limitation. Also, 
methods have been developed to overcome this (Mehta, 1978). A more 
critical flaw in the initial field is the approximately zero Reynolds 
stress uw . This stress controls both the transfer of turbulent energy 
from the mean flow, and the diffusion rate of the axial mean velocity. 
Clearly it is a critical quantity. Similarly, third order moments, 

■ - 2 ~ ■ 

e.g., u w , are also initially approximately zero. These control the 
self-diffusion of the turbulence. Another more subtle inadequacy in the 
initialization is the uniformity of the skewnesses and kurtoses . In a 
Tralee flow, the values for these quantities are usually larger near the 
wake edge due to the intermittancy of the flow there. However, the 
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uniform application of the form function J does not permit this 
intermittancy. A number of other inadequacies can also he pointed out. 

We have understood the nature of the limitations in the initial 
conditions. We concluded, however, that it would be useful to perform 
simulations with this scheme because (i) it would determine how signifi- 
cant these limitations were; (ii) it might show us how to better initial- 
ize the flow; and Ciii) it allowed us to study a flow which ultimately 
developed in approximately the proper manner. The results to be pre- 
sented below were obtained using this scheme from Orszag and Pao. We 
also experimented with several alternative initialization procedures. 
These alternative methods will he discussed in Section 4. 

In order to implement this scheme on the computer we had to finally 
choose a mesh size, a time step and a viscosity. We chose to treat a 
cubic box of side length L such that 


£ 


~ 8.0 , 


(61) 


This allowed for considerable spread of the wake before boundary effects 
became important. In x^ave number space, this resulted in a peak In the 
three-dimensional spectrum at 


st = 

and a cutoff wave number of 


iCp ” ^*^5 


(62) 


!<c $ ~ 


(63) 


which, was sufficiently large so that numerical truncation effects were 
not significant. The time step size At was chosen to he 



*" „ 07L & 


(.64) 
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and the viscosity V was selected so that the initial turbulent 
Reynolds number was 50, i.e. , 

— = so. css: 

y 

This is about the largest value that can accurately be treated on a 
32 x 32 x 33 mesh. 
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3. Numerical Results 

In this section we give the results of our direct numerical simula- 
tions of the axisymmetric, turbulent wake of a towed body. As mentioned 
in the previous section, we use the data of Pao and Lin (1973) at 
x/D - 36 , which is in the self-similar decay region, to initialize bur 
calculations. Tie present the results of the simulations of two independ- 
ent realizations, generated using different random number seeds in the 
initialization process. The combined results show both the capabilities 
of the numerical simulations to treat turbulent wakes, as well as provide 
an estimate of the statistical scatter in the results. We show first 
some constant contour plots of the instantaneous velocity and vorticity 
fields to help give a qualitative picture of the wake evolution. We 
then present plots of the decay of characteristic velocity and length 
scales, as well as total energies, with time. Finally, we present 
radial plots of various statistical quantities, and establish the degree 
to which the self-similar decay is modeled. 

3.1 Contour Plots 

Although statistical results are necessary to determine the capabil- 
ities of the simulations, it is useful to examine contour plots to see 
qualitatively the development of various features in the flow field. 
Figures 8 and 9 are contour plots of the longitudinal velocity fluctua- 
tion (u) and the magnitude of the fluctuating vorticity (to =y'G7oL ) 
respectively, in a typical plane perpendicular to the mean motion 
(x or x^) , at time zero. The vorticity plot especially shows the 
localization of the turbulence in the cylindrical region in the center.. . 
Note the axisymmetry of the vorticity field, as well as the manner in 
which it goes rather smoothly to zero with radial distance. These are 
both a result of the use of the smooth form function in the initializa- 
tion. Figures 10 and 11 are contour plots of u and to in a plane 
parallel to and passing through the x^ axis, at time zero. Again note 
the lack of bulges and other nonunif amities at the wake edge, especially 
in the vorticity field. Also, there appears to he no directional prefer- 
ence in the contours. This reflects the fact that the Reynolds stresses 
are initially approximately zero. 
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Figures 12 and 13 are similar contour plots of u and (jJ in 

the plane perpendicular to the x„ axis at nondimensional time 

1 A 

T (= of 7.73 , which is one- third of the way through the calculation 

The wake has grown somewhat, and the boundary has become much more 

convoluted. Also, some of the smaller scale features have been lost to 

molecular viscosity. Figures 14 and 15 are contour plots of u and 

to in the x^x,, plane through the x^ axis at T — 7.73 . Again note 

the convoluted nature of the wake boundary, and the growth of the wake. 

It is also apparent that the velocity field has developed a directional 

orientation. In the upper-half-plane, the constant contours tend to be 

oriented with their major axes in the second and fourth quadrants, 

whereas in the lower half-plane the contours are oriented with their 

major axes in the first and third quadrants. This can be shown to 

correspond to a positive Reynolds stress in the upper-half-plane and a 

negative stress in the lower-half-plane. Note that the mean velocity is 

positive in the direction of increasing x^ . The mean velocity 

gradient has kinematically distorted the Initially Isotropic-like field 

into one which has directional preferences . 

Figures 16 and 17 are contour plots of u and w in the x^x^ 
plane at time T = 15.46 , two-thirds of the way through our calculation. 
The continual growth of the wake, distortion of the wake boundary, and 
loss of small scale can again be observed. Figures IS and 19 give 
similar plots in the plane for this time. Again note the Increas- 

ingly convoluted surface, and the directional preference in the contour 
lines. The extent of convolution of the boundary regions of the wake 
will show up as increased values of the velocity (and velocity deriva- 
tive) skewness and kurtosis in those regions. 

3.2 Temporal Behavior 

We next turn our attention to the temporal behavior of characteris- 
tic velocity and length scales, and of the total energy. In Section 2.3 
we found that, in the self-preserving region, characteristic velocities 


We define the (constant) velocity % to be at time zero, and 

we define tR* to be r at time zero. 


m 


m 
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-2/3 


should decay as t w , characteristic lengths should grow as 
and both the total mean and turbulent kinetic energies decay as t 
In particular this implies that 


t 1/3 

L 5 

-2/3 
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(66) 

(67) 

( 68 ) 

(69) 

(70) 

(71) 


Here, U is the maximum value of the mean velocity deficit, u is 


m 


m 


the maximum value of the longitudinal tms velocity, r is the half- 


radius of the mean velocity, 
component of the turbulent intensity 


m 

i' T is the half -radius of the longitudinal 


E is the total (integrated) mean 


flow kinetic energy, E q is E at time zero, e is the total turbulent 
kinetic energy, t is a virtual origin in time, and C T 


T 


E 


U u m 

and C are constants. Ne will look for corresponding 


behavior in our simulations . 

Figure 20 is a plot of the temporal behavior of U_ and u„ for 


the two different realisations 


Consider U 


m 


m ■ ; ■ m 

first. Initially, the 


decay rate of is less than predicted by theory, or observed in the 

data. The dynamic equation for U is equation (32) : 


jt U + 7 fr(’’ uwr ') = r • 


( 32 ) 


*In presenting the simulation results, we show the horizontal axis in 
terms of both x/D (as in Pao and Lin), and also the nondimens ional 


% 


time T - , T = 0 corresponds to x/D = 36, and 

corresponds to x/D = 0 . 


T = T = -3.76 
o 
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The second term on the LHS represents the radial diffusion of mean 
momentum due to the Reynolds stress uw . The term on the RHS is the 
molecular diffusion term, and is small in this case. Initially, the 
mean velocity peak does not decay since uw = 0 . The Reynolds stress 
builds up rapidly, however, and by x/D - 90 1J is approximately 


.-2/3 


m 


following the t behavior for both realizations. Note that there 

is some difference in the results for the two realizations. 

~2 


The behavior of u is similar. 


m 


The equation satisfied by — 


is 


equation (34) : — 

^ Tp* — Air , i ,j L \ ~ - J, a&b 

• ~ v ™ m 


1,1 F z }) " * r i ; 


(34) 


The dynamics of u* are governed by a complex balance between energy 




feed from the mean flow (uw ■gpl 


the diffusion of u T by turbulence 


1 9 , u 2 . 

— -p— (r — w) 
r 3r 2 


the redistribution of u T between w r and v' 


, 1 ^P^ 

( - ? u 


and the viscous dissipation of u 1 .. Initially the generation terms and 
diffusion terms are both zero, due to the initialization scheme. Also, 
the dissipation is weak due to the initial lack of wave number amplitude 
correlations. (Hence, the wave number energy transfer is small, which 
is reflected, for example, in low values of the velocity derivative 
skewness.) However, by x/D of about 110, these effects all appear to 
have built up, and decays at roughly the proper rate for the two 

realizations. Some idea of the relative error can again be seen by 
comparing the results for the two realizations. 

Figure 21 is a plot of the temporal behavior of r and of r 


Considering first r 


m 


m *T 

, we . see that it displays the same effect as does 


IT . Initially it does not grow, but apparently after the Reynolds 

it increases at approximately the proper rate. The 
. Initially there is no growth, but by x/D equal, 
to about 110 it is growing at approximately the right rate in correspond- 
ence with the decay of u 

■ ■ xu 


stress has built up. 
same holds for r n 
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Figure 23 contains plots of E and e as functions of time. Multiply- 
ing equation (32) by pU and integrating over the crossplane gives the 
equation for E , i.e. , 


d 


.00 


= fr m 0 

o J 


o> 


(72) 


f/r -*frj>\) 
o 

The first term on the RHS represents the loss of mean kinetic energy to 

turbulence, while the second term is the molecular dissipation term, which 

is small. Mean flow energy is lost to the turbulence through the 
action of the Reynolds stress uw on the mean velocity gradient. The 
decay of E has established close to its proper behavior by about 

x/D = 100 . The equation for e is obtained by integrating the sum of 


2 2 

equation (38) and its counterparts for w and v 


and is 


i s = 




r uuj 


0 

w 


dr 


(733 


where £ is the viscous dissipation. Thus the time rate of change of 
e is a balance between energy feed from mean flow and viscous dissipa- 
tion. Again, the proper decay rate is approximately obtained by x/D 
about 110. The slightly slower decay rate could be due to wall effects, 
statistical scatter, the lack of sub-grid scale model, or an erroneous 
value for the initial integral scale A . 

Although the approximate proper decay rate has been attained by 
these quantities, their ratios have changed from the initial conditions. 
From Pao and Lin (1973) at x/D - 36 , 


Mm 


U = 31 ' 




jV 


2.0 


C74) 


and these ratios are maintained throughout the self -similar region. In 
our simulations, however, although we started the values of the above 
ratios properly (to within statistical scatter), the values in the self- 
similar range became 
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These values are somewhat different from the data, and are another 
reflection of our inadequate method of initialization. These ratios do 
vary, however, for different shaped obstacles. Our wake may be more 
representative of the walce of a sphere than a slender body. 

To summarize, initially the time development of the flow does not 
correspond to the self -similarity theory and the data due to inadequate 
initialisation, especially of the Reynolds stress. After the flow has 
adjusted, however, and the Reynolds stresses have built up, these 
quantities develop in time at approximately the proper self-similar 
rate, over a distance of about 100 to 270 diameters (os' a time of 

U m r T 

23.2 ^ . In the meantime the ratios — — and — change from their 

m m 

initial values , so that our calculations may be more representative of 
the wake of a different body than that used by Pao and Lin. 

3.3 Spatial Behavior 

In this section we present the radial distribution of various 

statistical quantities for different decay times. Self-similarity 

theory predicts that if the velocities are normalized by U (or u ) , 

m m 

and the lengths are scaled by r (or r ) , then the radial distribu- 

^ 1 r 

tion of a particular statistical quantity, say — Uf-y-) , should collapse 

m m 

for different times. The ability of this scaling to collapse the data, 
and the degree to which the resulting profiles agree with corresponding 
profiles from laboratory data is a good check of the capability of the 
numerical simulations to correctly model the physical processes occurring 
in the wake. 
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To. obtain the appropriate velocity and length scales, we assume 
that the flow is self-similar, and subjectively estimate best fits to 
the data given in Figures .20 and 21. The results are 



(The notation ( )* indicates a best estimate of the self-similar decay 
law for the numerical data.) These expressions are given by the solid 
curves in Figures 20 and 21. The following data are scaled using these 
expressions to compute the appropriate scaling quantity at the particular 
x/D . 

Figure 23 gives data from the first realization (Rl) for ^ 

m 

versus — for the locations x/D = 150, 205, and 260 ; U decreases 
r** m 

m 

by about 30 percent and r ■ increases by 20 percent over this range. 

We see, first of all, that the scaling collapses the data very .well, 

except near the origin, where the statistical scatter is expected to be 

the greatest. The laboratory results tend to follow the self-similar 

result rather well tFigure 6) , although the scatter is fairly large for 

x 

— > 1.5 , so that it is difficult to conclude that the agreement is 
m 

good in this region. Previous measurements of the wakes of two-dimensional 

bodies (see, e.g., Townsend, 1976) have consistently shown that the 

theory predicts somewhat larger values than are observed for — > 1.5 . 

m 

This is usually attributed to decreased .turbulent diffusion near the 

wake edge, due to the intermittent nature of the flow there. We see 

x 

that the simulation results tend to be slightly low for small z . This 

. ■■■ ■■ ■ m 



Flow Research Report No. 135 
March, 1979 

-32- 

is due to the fact that both realizations were talcen into account in the 

normalization. Data from the second realization lie slightly above the 

theoretical curve in this same region (Figure 24) . The simulation 

results also lie below the theoretical curve for — > 1.5 . The fact 

r** 

m 

that U initially started as an exact fit to the similarity result 
may indicate that the numerical simulations have added a realistic 
feature . 

Figure 24 is a plot of the same quantity at the same locations for 

the second realization (R2) . Again, the collapse is fairly good. The 

high values for small — indicates the data scatter to be expected 

m 

from one realization to another. The results for — ^ >0.7 are in good 

m 

agreement with those from R1 . We conclude from these results that the 
mean velocity behaves in a manner very close to the laboratory data, but 
that there is some scatter in the data near the origin. 

11 1 3T 

Figure 25 is a plot of — - versus — . for the downstream locations 

m m 

2$4d = 150, 205, and 260 for Rl . The dashed line is a rough fit of 

the laboratory data of Pao and Lin (1973). The laboratory data curve 

r T 

has been compressed so that - — for the data takes on the same value 

m 

given by equation (.75) instead of (74) i.e., 1.66. By plotting all the 

data versus — ^ , xje can compare the simulation results with laboratory 
m 

data, and simultaneously compare u 1 with U and other variables. The 
collapse of the simulation data is again fairly good. The scatter in 
the data for u r is expected to be significantly larger than that for 
U because the sampling error increases monotonically with the order of 
the statistical quantity (see, e.g., Kenney, 1954). For the choice of 
u* and r^ given by equations (77) and (78) the agreement with the 
laboratory data is only fair. If this data were matched independent of 
R2 , a better fit could be obtained. 

Figure 26 is a similar plot for R2 . Again the collapse of the 
data is fairly good except near the origin. This realization agrees 
much better with the laboratory data than the first realization. This 
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xs partly due to the choxce of u* and r* . When the data for — r 
r J m m u» 

m 

are superimposed for the two realizations, distinct differences in the 
two sets can he noted. This gives some idea of the magnitude in the 
scatter in the data to he expected in the simulation. 


Figure 27 is a plot of the normalized Reynolds stress 


uw 

u* 2 

m 


versus 


for 3</D = 150, 205, and 260 in 

m 


Rl 


The solid curve gives the 

prediction of the similarity theory, equation (31) , and the hatched 
lines gives the data of Pao and Lin (1973). The collapse is fairly good 
considering the error expected in computing the Reynolds stress. Further- 
more, the agreement with the laboratory data is reasonable. In Figure 28 
we present similar data from R2 . In this realization, the collapse is 
not quite as good. Also, a possible trend appears to be in the data, 
with the peak value increasing with x/D . Superimposing the two sets 
of data for the two realizations, shows reasonable agreement between the 
two. 


Figure 29 is a plot of 


2 

u w 

u * 3 

m 


versus 


r* 

. m 


for x/D = 150, 205, and 


260 for Rl 


2 2 

The radial flux of u is represented by u w and is an 


essential part of the dynamics of u (see equation (38)). The collapse 
is reasonable, except near the origin where the scatter in the data is 
expected to be greatest. Although we presently do not know of data with 
which to compare the result, we see that it has at least the proper 

~~2 

qualitative behavior; u is being diffused by the turbulence away from 


its peak (which occurs at about 


r -x 


m 


1 , Figure 25 ) , towards the origin 
r 


for — < 0.9 v and away from, the origin for — ^ > 0.9 . Figure 30 
m ~ m ~ 

gives the same data, for R2 . Here the -collapse is not nearly as good. 


Apparently this quantity has not fully adjusted in this range. Some of 
the scatter in the results could also be due to poor statistics. 

The intermittancy of the wake turbulence is reflected in the flat- 
ness factor of the velocity field. Townsend (1949) has shorn for wakes, 
and Wygnanski and Fiedler (1970) for mixing layers, that the velocity 
kurtosis is approximately equal to the inverse of the intermittancy 
factor. Figure 31 is a plot of the radial distribution of the velocity 
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Icurtosis for R1 and the same three dox-mstream locations. At the start 
of the calculation, the profile is flat (Table 1), due to the method of 
initialisation, indicating there is no intermitt ancy in the flow. As 
the walce develops, however, the outer regions of the wake become highly 
contorted (e.g. , Figures 8 to 19), and the velocity becomes highly 
intermittant, which is reflected in the sharp increase in the leurtosis 
at the wake edge. The dashed line is a rough fit of the data of Pao and 
Lin (1973), and the agreement between the simulations and laboratory 
data is reasonable. 

Figure 32 is the same plot for R2 . While the collapse of the 
data is good, the agreement with the laboratory data and with Rl is 
poor, especially in the locations of the peak value. Apparently in R2 
the walce edge is neither as contorted nor as far out as in Rl . This 
again is an indication of the deviations to expect from one realization 
to another. To compute the kurtosis accurately, xxre may need to compute 
seyeral realizations and use ensemble averaging. 

The velocity skewness is also an indicator of the intermittancy of 
the flow field. Figure 33 is a plot for the velocity skewness for Rl 
for the same three downstream locations. Again the collapse is fairly 
good, and the wake edge is clearly distinguishable through the large 
skexmess values. The agreement with the laboratory data is reasonably 
good. Figure 34 gives the same plot for R2 . The results are not in 
good agreement with either the laboratory data or with the skewness data 
from Rl . This is another indication that the wake edges for the two 
realizations are in different radial locations. 

To summarize the radial plots, in general the collapse of the data 
using the self -similarity scaling is very good. Also the agreement with 
laboratory data is generally good, except for two cases. The first is 

the turbulence intensity (u r ) profile, where the ratio — has changed 

m 

from its initial value, equations (74) and (75) . However, when u T is 

r 

r T 

plotted against — , or, equivalently, when the ratio of — - given by 
r T r m 

(75) is used to plot the laboratory data, the agreement is then fairly 

good. Secondly, some significant differences were observed in the 

skexmess and kurtosis for the second realization. 
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4 , Other Methods of Initialization 

It is obvious from the simulation results presented that our method 
of initialization, although properly treating many important turbulence 
properties, was inadequate to properly initialize all of the significant 
turbulence quantities in the wake. This inadequacy was particularly 
exhibited in the initial Reynolds stresses, higher order flux terms, and 
skewnesses and kurtoses. ’ 

We have begun a preliminary exploration of ways to improve our 
initialization procedure. Although we have not yet had time to reach 
definite conclusions on better ways to initialize the turbulence in the 
wake, we thought it would still be useful to discuss the results of 
these efforts. For the successful utilisation of direct numerical 
simulations of turbulence, the difficulties inherent in properly ini- 
tializing turbulent shear flows may be as important to resolve as the 
questions arising from subgrid-scale modeling. 

We think that the key improvements to be made are in the proper 
initialization of the Reynolds stresses and higher order flux terms, 
such as the turbulent flux of turbulent energy. As mentioned in the 
last section, these terms represent the diffusion of the mean velocity 
by the turbulence, the generation of turbulent energy, and the self- 
diffusion of the turbulence. Improper initialization eliminates the 
possibility of an exact initialisation of the wake. A secondary improve- 
ment of importance is the correct specification of the structure of the 
flow at the wake edge, which is reflected in the skewnesses and kurtoses. 

A possible method to treat these problems is the extension of a 
technique used by Rogallo (1977) and Mehta (1978) in computing homoge- 
neous turbulence decay.. As our calculations proceed in time, the Rey- 
nolds stresses and higher moments do build up to approximately their 
proper values. During this buildup period, however, the mean velocity 
has diffused outward and the turbulence has decayed somewhat, so that 
all of the proper conditions are not attained. The idea is to reini- 
tialize the velocity field by retaining the phase relationship in each 
Fourier mode that has been attained when the stresses have built up, but 
reimposing the initial Fourier amplitude. Since the information about. . . 
the Reynolds stresses and higher movements is mainly contained in the 
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phase relationships of the inodes, it was hoped that the initial Reynolds 
stress behavior would be desired. Because the energy levels were read- 
justed to their initial values, it was also hoped that the initial 
profiles of turbulence intensity would again take on their proper values. 

We implemented this technique on the computer. Two free parameters 
at our disposal were the number of time steps between reimposing the 
initial Fourier amplitudes (N^) , and the number of times we reimpose 
the initial amplitudes (N) . It takes about 100 time steps for the 
Reynolds stress to attain its proper value. We experimented with the 
following values of and N : (H, - (1» 120), (I* 80), (2, 40), 

(5, 20), (10, 10). In general, we found that the best overall results 
were obtained for cases with smaller N . We will present the results 
for (W, N ) - (2, 40) which is representative of the better results, 
but which also points out some of the difficulties with the method . 

Figure 35 contains a plot of the initial profile of the longi- 
tudinal turbulent intensity u* before the field was reinitialized. 

Figure 36 contains a plot of the corresponding Reynolds stress uw . 

The turbulent field for this realization was reinitialized using values 
of N = 2 and W T = 40. The reinitialized quantities are also shown in 
Figures 35 and 36. The new Reynolds stress has developed to the proper 
form, although its peak value is somewhat larger than the corresponding 
value in the laboratory data. The difference is typical of the statisti- 
cal fluctuations that we have seen in uw in our simulations (Figures 27 
and 28) . The new turbulence intensity profile has changed substantially 
from the initial profile. Some of the energy near the axis has moved to 
outside the wake. Apparently, although the proper phase information for 
the Reynolds stress has been retained, the phase information necessary 
to maintain the proper profile for u’ has been degraded. Because of 
this behavior, we find this method as presently used unacceptable, and 
are working on its improvement. 

We have also experimented with several other methods. The first 

is a shooting method, which takes into account the fact that, although 

the conditions ultimately attained are not exactly as desired (e.g., 
u 

jj— = 0.47 instead of 0.32), they are close. Thus it may be possible to 
m 

adjust the initial field to shoot to the proper values for these quantities 
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We have implemented this shooting method and are presently testing it. 
Another method* suggested by Main (.1978) , involves the use of the eigen- 
functions of the most unstable modes for the wake stability problem. The 
idea is that the flow structure which supports the Reynolds stresses may 
be similar to these unstable modes. Their Reynolds stress can be easily 
computed, and the eigenfunctions added into our initial field to produce 
the proper stresses. We have yet to implement this method. 

Several other ideas are under consideration. We are still in the 
process of studying the applicability of these various methods, with the 
ultimate goal of the proper initialization of the turbulence field. 



Flow Research Report No. 135 
March, 1979 


-38- 

5, Conclusions 

From the numerical simulation data presented in the preceeding 
section, we can draw a number of conclusions about the capabilities of 
the direct numerical simulations of turbulence. In general, the agree- 
ment of the numerical simulation results with the laboratory data and 
the self -similarity theory was good. In terms of the temporal behavior 
of the simulated flow over a span of time of about 15 ^ , 

(i) the maximum mean and turbulent velocities decayed at 

—2/3 

approximately the proper rate ((x/D ) ; 

(ii) both the mean and turbulent walce radii increased at 

-i/3 

approximately the proper rate ((x/D) ) ; 


(iii) 


the total integrated mean and turbulent energies decayed 

-2/3 

at approximately the proper rate (Ox/D ). 


Furthermore, over this same interval of time, using self -similarity 
scaling, we were able to obtain reasonably good collapse, and also 
reasonable agreement xvith laboratory data, for 


(i) the mean velocity profile 

(ii) the turbulent intensity profile 

(iii) the Reynolds stress profile 

~ 2 ~ 

(iv) the profile for u w (no laboratory data was available) , 
and 

(v) the profiles for skewness and kurtosis. 
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In addition, several realistic features, which were not in our initial 
conditions, subsequently appeared later in time In the numerical simula- 
tions. These include 

(i) the development and maintenance of the Reynolds stress, 

Cxi) the development and maintenance of higher order correla- 

2 

tions, e.g. , u w ; 

(Hi) the development of intermit tency near the wake edge, 

which qualitatively shows up as bulges In the walce, and 
is reflected by the development of the skewness and 
kurtosis to reasonable values. 

It also appeared that a subgrid-scale model is not a primary necessity 
for computing this free shear flow, although a subgrid-scale model may 
prove to he necessary in the future when more careful comparisons with 
data, especially for the decay of turbulent energy, are made. 

The comparisons were hampered by two inadequacies in the present 
calculations. The first was the initialisation process. The turbulent 
field, especially in terms of the Reynolds stresses, was not adequately 
initialized, which limited the conclusions that we could draw concerning 
numerical simulations. Secondly, there was a significant amount of 
scatter in our data, both, in an individual realization, and from one 
realization to another. 

Given these limitations, the present work greatly increases our 
confidence in the capabilities of direct numerical simulations to 
accurately treat the physics of turbulent shear flows. The turbulence 
in the axl symmetric wake is governed by a complex balance among the 
following mechanisms; generation of energy from the mean flow, diffu- 
sion of turbulence radially by the turbulence, inter component transfer 
of energy, and dissipation. All of these processes have to be repre- 
sented faithfully for the resulting statistics to possess the proper 
behavior. (It should be noted that there are no adjustable constants in 
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our simulations.) Thus we feel that the outlook is good for using direct 
numerical simulations to predict statistical results which depend on the 
large-scale features of free turbulent shear flows. 

Prom our work, the emphasis for future research work in this field 
has become more clear. The following are topics which we think should 
be addressed in the near future: 

(i) The method of initialization should be improved to 
include the proper initialisation of the Reynold stress 
and higher order moments. Until this is done close 
comparisons with data will be difficult. (The initializa- 
tion of the turbulent velocity field for a statistically 
time dependent problem is analogous to the specification 
of turbulent inflow boundary conditions for spatially 
nonhomo geneous problems.) 

(ii) The present simulations should be examined in more depth 
to look, for example, at the behavior of spectra and 
correlations. At the same time, the effects of a subgrid- 
scale model should be examined. 

(iii) Because of the limitations in statistics obtained from 
these 32 x 32 x 33 grid calculations, it is probably 
necessary to use a 64 x 64 x 65 grid and/or compute a 
number of realizations in order to obtain statistics 
which will adequately test the method. 

(iv) It is important to develop techniques for specifying 
proper inflow and outflow boundary conditions to enable 
the treatment of more general problems . 
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Table 1, Initial Data - Realization 1. 


n 

Cylindrical 
. Annulus 

N 

n 

Number of 
Points 


r/r^ 

. m 

Cr - 
m 

0*82815) 

U 

u' 

(axial) 

w’ 

(radial) 

V* 

(angular) 

2 

u w 

(X10 -2 ) 

3 

w 

(X10 -2 ) 

2 

V v 

(X10 -2 ) 

1 

32 

0.1768 

0.1185 

1.0000 

0.3588 

0. . 

0. 

0. 

0. 

0. 

2 

256 

0.0625 

0.3556 

0.9363 

0.3401- 

0.3714 

0.4055 

0.3207 

0.6609 

0.6170 

3 

512 

0.0442 

0.5927 

D. 7866 

0.3370 

0.3688 

0.4054 

0.7777 

-0.4200 

0.1491 

4 

640 

0.0395 

0.8298 

0.6178 

0.3254 

0.3624 

0.3619 

-0.1412 

-0.4814 

0.5297 

5 

76 B 

0.0361 

1.0669 

0.4537 

0.3037 

0.3525 

0.3308 

0.0781 

-0.0756 

0.4250 

6 

1280 

0.0280 

1.3040 

0.2884 

0.2842 

0.3113 

0.2956 

0.1706 

0.0645 

0.0712 

7 

1152 

0.0295 

1.5411 

0.1717 

0.2748 

0.2636 

0.2579 

0.0835 

0.1191 

-0.0121 

8 

1536 

0.0255 

1.7782 

0.0947 

0.2046 

0.2146 

0.2086 

0.0322 

0.0426 

0.0301 

9 

1792 

0.0236 

2.0153 

0.0462 

0.1683 

0.1656 

0.1719 

0.0124 

-0.0227 

0.0004 

10 

1792 

0.0236 

2.2524 

0.0212 

0.1402 

0.1279 , 

0.1386 

-0.0135 

-0.0226 

0.0145 

11' 

2176 

0.0214 

2.4895 

0.0090 

0.1035 

0.0963 

0.1030 

-0.0020 

-0.0104 

0.0119 

12 

■2048 

0.0221 

2.7266 

0.0036 

0.0761 

0.0728 

0.0783 

: -0.0025 

0.0042 

0.0036 

13 

2560 

0.0190 

2.9637 

0.0013 

0.0528 

0.0529 

0.0612 

-O.OD14 

0.0028 

-0.0006 

14 

2944 

0.0184 

3.2008 

0.0004 

0.0369 

0.0369 

0.0450 

-0.0004 

0.0003 

-0.0005 

15 

2B16 

0.D188 

3.4379 

0.0001 

0.0285 

0.0260 

0.0352 

-0.00D1 

0.0000 

0.0000 

16 

3072 

0.0180 

3.6750 

0.0000 

0.0218 

■ 

0.0175 

0.0252 

0.-0000 

0.0000 

0*0000 












Table 1, (Cont.) 
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Figure 8. Velocity Field Contour Plot 
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Figure 12. Velocity Field Contour Plot 
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Figure 13. Vorticity Magnitude Contour Plot. 
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Figure 15. Vorticity Magnitude Contour Plot 
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Figure 16. Velocity Field Contour Plot 
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Figure 17. Vorticity Magnitude Contour Plot 
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Figure 18. Velocity Field Contour Plot. 
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Figure 20. Decay of Maximum Mean Velocity (U m ) and Maximum Axial Turbulent Intensity (u m ). 
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Figure 21. Growth of Mean Wake Radius (r m ) and Turbulent Wake Radius (r T ). 
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Figure 22. Decay of Integrated Mean and Turbulent Energies. 
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re 23. Normalized Profile of Axial Mean Velocity for Realization 1. 
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Figure 25. Normalized Profile of Axial Turbulent Intensity for Realization 1. 
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Figure 26. Normalized Profile of Axial Turbulent Intensity for Realization 2. 
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Figure 27. Normalized Profile of Reynolds Stress for Realization 1. 
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Figure 28. Normalized Profile of Reynolds Stress for Realization 2. 
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Figure 29. Normalized Profile of Radial Flux of Axial Turbulent Energy for Realization 1. 
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Figure 30. Normalized Profile of Radial Flux of Axial Turbulent Energy for Realization 2. 
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Figure 31 . Profile of Velocity Kurtosis for Realization 1 . 
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Figure 32. Profile of Velocity Kurtosis for Realization 2. 
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Figure 33. Profile of Velocity Skewness for Realization 1 . 
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Figure 34. 


Profile of Velocity Skewness for Realization 2. 
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Figure 35. Initial Axial Turbulence Intensity. 
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